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Abstract 

Wc describe a systematic method for complete enumeration of configuration 
classes (CCs) of the spin-1/2 Ising model in the energy-magnetization plane. 
This technique is applied to the antiferromagnetic Ising model in an exter- 
nal magnetic field on the square lattice, which is simulated using the tomo- 
graphic entropic sampling algorithm. We estimate the number of configurations, 
n(E,m, L), and related microcanonical averages, for all allowed energies E and 
magnetizations m for L = 10 to 30, with AL = 2. With prior knowledge of 
the CCs, we can be sure that all allowed classes are sampled during the simula- 
tion. Complete enumeration of CCs also enables us to use the final estimate of 
Q{E,m, L) to obtain good initial estimates, flQ{E,m, L'), for successive system 
sizes (L' > L) through a two-dimensional interpolation. Using these results 
wc calculate canonical averages of the thermodynamic quantities of interest as 
continuous functions of temperature T and external field h. In addition, we 
determine the critical line in the h-T plane using finite-size scaling analysis, and 
compare these results with several approximate theoretical expressions. 

Keywords: Ising model; antifcrromagnct; configuration classes; Monte 
Carlo simulation; phase diagram. 
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1 Introduction 



The most important task of cntropic sampling algorithms [I] ~ [6] is to visit the 
full configuration space (CS) to obtain good estimates of the number of config- 
urations, n, as functions of the energy E and other quantities of interest. In 
studies of the Ising model in an external field, for example, we require m) 
with m the magnetization; each allowed (E, m) pair defines a class of configura- 
tions (CC). Only if we know beforehand the possible values of {E, m) for a given 
system size, can we be sure that all CCs are sampled during the simulation. 

Figure [T] shows the CCs in the n — m plane for the spin-1/2 Ising model 
with nearest-neighbor (NN) interactions, on a square lattice oi L x L sites with 
periodic boundaries. [We use n to denote the number of NN pairs of spins with 
the same orientation; the interaction energy of the antifcrromagnctic (AF) Ising 
model is E — — 2(L^ — n).] Although the CCs tend to fill in a triangular region, 
some gaps are evident near the lower apex and along the upper edge. Knowing 
just which (n, m) values arc allowed for a given lattice size is important if we 
are to implement cntropic sampling with confidence. In this paper we present 
a method for systematically enumerating all CCs of Ising models in the (n, m) 
plane. 



2 -.^ .■ - 

1.5 - 'lllllllllllllllllllllllllllllllllll^ 

0.5 '■ 

[ ^ ^ ^ ^ 

-1 -0.5 0.5 1 

V 

Figure 1: Allowed configuration classes for system size L = 10; rj = n/L'^ and u = m/L'^. 

We study the spin-1/2 antiferromagnetic Ising model in an external magnetic 
field, whose energy is given by 

N 

H = — J (JiCTj — /i^^ cTi = —E — hm, (1) 

<i,j> 1=1 

where ct^ = ±1, < i, j > indicates a sum over NN pairs of spins, h is the external 
field, and N is the number of spins; the model is defined on a square lattice of 
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L X L sites, with periodic boundary conditions. Unlike the ferromagnetic Ising 
model (J > 0), which exhibits a unique critical point in the h — T plane and 
has an exact solution [7], the AF model (J < 0) possesses a critical line, which 
is not completely understood. 

Various approximate methods have been applied to determine the critical 
line of the AF Ising model on the square lattice IH]~IIH]; these results, however, 
do not agree altogether. Binder and Landau [19] estimated the critical line via 
Monte Carlo simulation, obtaining very good agreement with the approximate 
closed-form expression of Miiller-Hartmann and Zittartz [5], raising the possi- 
bility that the latter expression was in fact exact. (It was later shown that this 
is not so [H].) Hwang et al. [ID] studied the AF Ising model on the square 
and triangular lattices using a microcanonical transfer matrix method and the 
Wang-Landau algorithm [3] . They performed an exact enumeration of the num- 
ber of configurations, VL{E,Tn), and found CSs in the E — m plane similar to 
that shown in Fig. [TJ 

Using the tomographic entropic sampling (TES) algorithm [B] we estimate 
fl{n,m, L), and associated microcanonical averages, for lattice sizes L = 10 to 
30, with an increment AL ~ 2. We then calculate the canonical averages of 
the thermodynamic quantities of interest. Using these data we map out the 
critical line in the h — T plane, and compare our results with several theoretical 
expressions. 

Prior determination of the set of allowed CCs is an important tool to verify 
the quality of the sampling: we want to be sure that all CCs are visited during 
the simulation. Since this algorithm uses an initial guess, ^{){n, m, L), to begin 
the study, it is convenient to use the final estimate r2Ar(n, m, L) (after the iV-th 
iteration) to obtain the initial guess il,o{n,m, Lq) of the next system size to be 
studied {Lq > L). As we will show, good initial estimates, nQ{n,m, L), can be 
obtained using a two-dimensional interpolation because we know a priori the 
set of allowed CCs. 

This paper is organized as follows. In Sec. [2] we define the basic CCs and the 
respective allowed values of (n,m) for the spin-1/2 Ising model on the square 
lattice; the main goal of this section is to find all gaps in the (n, m) plane. In 
Sec. m this information is used in simulations of the AF Ising model via the 
TES algorithm. There we describe the method used to determine ^^{n, m, L2) 
via two-dimensional interpolation of the final estimate, r2Ar(?T., m, Li), of the 
previous system size studied. Simulation results are reported in Sec. HI for the 
order parameter and the staggered susceptibility as functions of h and T. Points 
along the critical line in the h — T plane arc obtained using finite size scaling 
analysis, and the results compared with several theoretical expressions. We 
summarize our findings in Sec. [5] 

2 Configuration Classes 

As pointed out above, one of the main problems in entropic sampling methods 
is the prior determination of the complete set of configuration classes for a given 
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system size. Let us denote by iV+ and iV_ the number of up and down spins, 
respectively, on a square lattice with N = + = spins. The number 
of pairs of opposite spins, it, and are related to n and m, respectively, via 

n = 2L^ -u (2) 

and 

TO = 2N+ - L^. (3) 

Thus, once the possible values of (iV+, u) are determined, so arc those of (ri, m). 
Note that u, n and m can only take even values. The gaps in CS fall near 
the maximum and minimum values of n {n„iax and rimin, respectively) for a 
given TO. Therefore, we will identify the possible values of u near its maximum, 
Umax, and minimum, u,nin , for a given N-^ <E [0, -^]; note that the number of 
configurations is symmetric under interchange of N-^- and . 

2.1 Determining cind nearby classes 

2.1.1 Compact configurations 

Compact configurations consist of a square or rectangular cluster with up 
spins. To begin, consider the case of a square cluster of size I x I {2 < I < L) 
with = l^, as is illustrated in Fig. [2l It is evident that this configuration 
corresponds to the minimum value of u, u^"^ = 41. A configuration with the 
same number of up spins, and u = Af+ + 2 is obtained by transferring an up spin 
from one of the corners of the square to an edge, as shown in Fig. |31 From this 
configuration, further rearrangements leading to u = Umin + 'i-, etc., are possible. 
When is not a square number, the most compact configuration (i.e., with the 
smallest perimeter) is a rectangle, or a square or rectangle with an incomplete 
layer of sites along one edge. For iV+ = — we have Umin = 4/ — 2, while for 
l{l — 1) < VV+ < l^, u„iin = 4L In all cases, moving a corner site to an edge, one 
constructs a configuration with u = Umin + 2, and further arrangements yield 
additional increases in u. 

2.1.2 Extended configurations 

Thus there arc no gaps in the largc-n region due to compact configurations. Such 
compact configurations, however, are not necessarily the ones that minimize u 
for a given value of N+. Consider for example the case = kL with k an 
integer. The cluster of up spins can be arranged to wrap around the lattice, 
yielding Umin = 2L, independent of k. We call such a configuration extended. 
For k greater than a certain value, on the order of i/4, the configurations that 
minimize u arc extended rather than compact. Suppose now that for a given 
N^, Umin is obtained with an extended configuration, and that all compact 
configurations have u strictly greater than Umin- If lies between kL and 
[k + 1)L, then the configuration that minimizes u has at least two corners, and 
by the previous construction, configurations with u = Umin + 2 exist. But if 
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Figure 2: Basic compact configuration. Up and down spins are represented by "x" and 
"•" , respectively; wavy lines represent pairs of opposite spins. 
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Figure 3: Modified compact configuration. The new up and down spins are represented 
by "(Xi" (previously "•") and "0" (previously "x"), respectively; double straight and double 
wavy lines represent new pairs of identical and opposite NNs, respectively. 
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iV+ = kL, the minimizing configuration has no corners and any modification 
yields a configuration with u > Umin+^- This is how the gaps near the maximum 
values of n arise. 

Summarizing, if L and are such that Umm is obtained with a compact 
configuration, then there are configurations with Umin + 2, Umin +4, . . . , etc., 
and no gap exists. If, on the other hand, Umm is realized only for an extended 
configuration, and iV+ is an integer multiple of L, then there arc no configura- 
tions with U = Umin + 2. 

2.2 Determining Umax and nearby classes 

The largest possible value of w, Umax, occurs in a configuration such that A'^-i- = 
LV2 with spins arranged in a chess board (CB) configuration, such that all 
up spins have down spins as NNs and vice versa. One readily verifies that for 

< N+ < i^/2, the maximum number of unlike NN pairs is Umax = 4A^+. 
Starting from the CB configuration, we can reduce u by exchanging an up and 
a down spin. If the exchanged spins arc NNs, the resulting configuration has 
u = Umax — 6; otherwise one has u ~ Umax — 8. Thus, for — L^/2, there 
are no configurations with u = Umax — 2 or Umax — 4. One readily verifies that 
configurations with u = Umax — 10, Umax — 12, etc., can be obtained via further 
exchanges of spins. 

For iV_|_ = — 1, we have Umax = 2L^ — 4. Such configurations can be 
constructed by flipping one up spin in the CB configuration. Starting from this 
configuration, exchanging a pair of spins, one can reduce m by 4, 6 or 8, but 
there is no rearrangement which reduces u by just two. Thus for ~ — 1, 
there is no configuration having u = Umax — 2; configurations with u = Umax — 4, 
Umax — 6, etc. do exist. 

Finally, we note that for 1 < A^-|- < L? /2 — \, there are no gaps in the 
neighborhood of Umax- To verify this, consider a configuration obtained by 
flipping k up spins in the CB configuration, so that = — fc, where 

1 < k < By hypothesis there are now at least four more sites with down 
spins than with up spins. Starting from a configuration with u = Umax, in which 
each up spin is surrounded by down spins, we can create a single NN pair of up 
spins, with all six neighbors down, and with all remaining up spins completely 
surrounded by down spins. In this manner, u is reduced by two. Configurations 
with u ~ Umin " 4, Umin ~ 6, ctc. Can bc obtained by further exchanges of up 
and down spins. 

Using the facts summarized above, it is straightforward to construct an 
algorithm that determines which values of u are possible, for a given L and N+, 
and thereby which values (ri, m) arc accessible for a given system size. In the 
simulations reported below, we have verified that our entropic sampling scheme 
converges to visit all allowed classes. 
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3 Implementation 



Using tomographic sampling, wc study the antifcrromagnctic Ising model in 
an external field on the square lattice; we consider periodic boundary condi- 
tions and NN interactions. The CCs of the systems are defined in the energy- 
magnetization space (n, m). 

The TES method is applied in order to generate estimates of il{n, m, L). For 
the smallest system size {L = 10) we begin with a guess, r2o(n, m), obtained 
using a mean-field approximation. For subsequent system sizes, however, we use 
a two-dimensional interpolation of UnIti, m, L) (the final estimate of (n, m) for 
the smaller system size, L) to construct no{n,m, L'), where L' > L. For most 
studies we use five iterations, each one with -/Vsim ~ 10 initial configurations, 
which are simulated for Njj = 10^ lattice updates or Monte Carlo steps. 

Let us denote by T{n, m) and r'(no = n + An, mp = m + Am) the CCs that 
contain configurations C and C, respectively. The simulation uses a single spin- 
flip dynamics, so that the possible variations of n and m are An = 0, ±2, ±4 
and Am = ±2. At iteration j, the acceptance probability for the transition 
(C ^ C) is 



Pj{T F') = min 



f],_i(r) 



(4) 



f},_i(F')' 

These probabilities are stored in a table. For each configuration generated, be 
it a new one (if it is accepted) or the same (if it is rejected), we update the sums 
used to calculate the microcanonical and canonical averages of ((>, |0|, 0^, and 
(p'^, where 

= niA - ms (5) 

is the order parameter (staggered magnetization); rriA.B sltc the magnetizations 
of the two sublattices. At the end of each iteration j the estimate of D,j{n,m) 
is refined according to 

^A^) = ^,^j-i{n (6) 

where Hj{T) is the histogram containing the number of times the CC F is visited 
during the sampling, and Hj{T) is the average of Hj{T) over all accessible 
CCs; the acceptance probability [Eq. (|4])] is updated using the new estimate of 
r2j(n, m) and the histogram, Hj{n,m), is set to zero. 

A single iteration of our method consists of ten independent simulations, 
each involving 10^ lattice updates, and each beginning from a different initial 
configuration. (By a lattice update we mean one attempted flip per spin; the 
initial conflgurations include, high, low, and intermediate interaction energies 
and both signs of the magnetization.) 
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3.1 Determining riol?^) mean field approach 



As noted above, for the smallest system size we use an estimate of rio(n,m,L) 
obtained via a mean-field approximation, specifically 



1 1 (v — (v))'^ 

where cr = a^m) = ■\/var(n), and 



In a 



const. 



(7) 



lnf2(m) 
L2 



In 2 - -[(f + m/L^) ln(l + m/L^) + (f - yn/L^) ln(l - m/L^)]. (8) 



To obtain this expression, we first note that D,(m, L) = '^ji^{n,m, L) = 
(_/v+)' ^^'^ Stirling's approximation. The dependence on n is then obtained 
by estimating (n) and var(n) for a given m and L, using a random- mixing 
approximation, and supposing that n follows a Gaussian distribution. 

In Fig. Uwe plot the estimate of \nilQ{n,'m,, L — 10) given by Eq. ([7]). We 
present in Fig. [S] the final estimated value of lnri5(ri, m, L — 10) after the 5th 
iteration of the simulation; this result is quite similar to that presented by 
Hwang et al. |20j for the square lattice, which was obtained using Wang-Landau 
sampling [3] . We note that the differences between the initial estimate obtained 
via mean field approximation and the final simulation result of lnr2(n,m,L = 
10) are more evident near the edges of the CS, specially near to the maximum 
values of n. To have a better idea of how close this initial estimate is to the final 
simulation result, we plot in Fig.[n]the difference between that estimate and the 
final result of the simulation for L = 10. Analyzing this figure, it is again clear 
that the differences are larger near the edges of the CS. 
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Figure 4: Estimate of In ^o{n, m, L = 10) via mean-field approximation. 
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Figure 5: Final estimate of lnf25(n, m, L = 10) after the 5i/i iteration of the simulation. 




Figure 6: Difference between lnf2o(n, m, L = 10), estimated via mean- field approximation, 
and the final simulation result after the 5th iteration, lnQr^{n,m, L = 10). 



3.2 Determining no{n,m, L): interpolation 



We expect that the closer the estimate ^o{n, m) is to the (unknown) exact value, 
fie(n, m), the faster the simulation will converge, and the more accurate our final 
estimate will be. Since the mean-field estimate worsens as the system size grows, 
we only use it for the smallest system size studied; initial estimates of subsequent 
system sizes are obtained by interpolating the final estimate, J7Ar(n, m, L), of 
the previous system size studied. Our procedure is based on the existence of the 
limiting microcanonical entropy density as a function of the intensive parameters 
77 and v (the bond and magnetization densities, respectively), 

s{ri,v) = lim lnri(nL, tol, L), (9) 

L— foo h 

where tll ~ i]L^ and nriL — vL? . (Since n and m arc restricted to even integers 
we have = r/L'^ + 0(1/L^) and similarly for mL.) The idea is then to write 

-i^ lnl^o(^',m',i') = -llnOjv(??L',i^i',i), (10) 

where r] = n' /L'^^ v = m'/L'^, and the r.h.s. is evaluated by extending IuJItv to 
noninteger n and m via extrapolation and interpolation. Using this approach, we 
obtain better estimates, as is shown in Fig.[71 (Note that the largest differences 
continue to faU along the edges of the CS.) 




Figure 7: Difference between \nQ,o(n, m, L = 18), estimated by interpolating tfie final result 
of L = f6, and the final simulation result after the 5th iteration, lnQr^{n,m, L = 18). 



3.3 Extrapolation and interpolation 

Suppose we wish to construct the initial estimate flo{n,m, L2) on the basis of 
the simulation results for a smaller system, r2Ar(n, m, Li). We could do this via 
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interpolation if every CC of the larger system were surrounded by four points 
of the smaller one in the 77-/X plane. Fig. |8] shows, however, that along the 
edges of the CS, the points corresponding to classes of the larger system are not 
surrounded by four points of the smaller one. For those points, one could in 
principle use extrapolation rather than interpolation. We found, however, that 
direct extrapolation yields poor estimates for fi. 
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Figure 8: Region of the (rj, v) plane with every possible configuration classes for system 
sizes L\ = 26 and L2 = 28. 

We obtain better estimates by first extrapolating the points along the edges 
of the CS for the smaller system, such that every point of the larger system is 
surrounded by four points of the smaller. Figure |9] shows the CS with accessible 
and extrapolated CCs for L = 10. Following this extrapolation we perform a 
linear two-dimensional interpolation as per Eq. ([T0| . 



4 Results 

In this section we present results of the AF Ising model on the square lattice 
in an external field; periodic boundary conditions are employed. We use TES 
to simulate systems of sizes L = 10 to 30, with AL = 2. To calculate the 
uncertainties we perform five independent studies for each system size. We 
plot in Fig. [TU] the staggered magnetization (order parameter) per site, ^, as a 
function of /i at T = 0.2. We can see that 4> decreases considerably between 
h = 3.85 and h ~ 3.90; this behavior suggests a critical point, he, marking a 
phase transition from the AF to the paramagnetic state. Figure [TT] shows the 
staggered susceptibility per site, 

^ ^{{<p') - m, (11) 



11 



as a function of /i at T = 0.02. As L grows the peaks tend to the critical point, 
he- The specific heat per site, c, (not shown) has a similar behavior. 



4.1 Phase diagram 

Using finite size scaling analysis [21] we estimate the critical line, hc{T), or, 
equivalently, Tc{h), via the relations: 

hc{Y,nax,T, L) = K{Y^ax.T) + a^jL + aa/i', (12) 

and 

T,{Yrnacc,K L) = re(r™,„ h) + / L + ^a/i', (13) 

where hc(Ymax^T^L) is the field at which Y (the specific heat or the staggered 
susceptibility) takes its maximum for a given temperature and system size; 
Tc{Ymax^ h, L) is defined in an analogous manner. The estimates obtained using 
the maximum of c and (f> are averaged to yield hc{T) and Tc{h); Figure [1^ illus- 
trates the procedure for estimating he, for T = 0.02. Our estimated points for 
the phase boundary arc shown in Fig. [T3] along with the approximate expression 
derived by Miiller-Hartmann and Zittartz [8]: 

c„si(i).,i„„'(|)^ (14) 

Our simulation results are in good agreement with their expression. For the 
critical field, the greatest relative difference between theory and simulation is 
about 0.9%, which occurs at Tc ~ 1.8. 

In Table [T] we compare our simulation estimates for the critical magnetic 
field hc{T) with some theoretical approximations. For temperature T < 1 we 
find good agreement with the estimates of Monroe [T7] (whose analysis involves 
a free parameter w), Wu and Wu [12], and Blote and Wu [13], whereas there are 
significant discrepancies in relation to the other approximations. At higher tem- 
peratures, differences appear between simulation and the predictions of Monroe, 
Wu and Wu, and Blote and Wu. These may reflect a small systematic error or 
an underestimate of simulation uncertainties. We intend to examine this issue 
in greater detail in future work. 



5 Conclusions 

The complete enumeration of CCs is of fundamental importance to study the 
antifcrromagnetic Ising model using tomographic entropic sampling. The deter- 
mination of CCs for entropic sampling of Ising models also enables us to obtain 
good initial estimates for the configuration numbers ^lo{n,m, L), using two- 
dimensional linear interpolation; the initial estimate is reasonably close to the 
final estimate ^n- Despite the relatively small system sizes used in this study, 
we obtain a good estimate for the critical line in the temperature-magnetic field 
plane. Further details on critical behavior will be published elsewhere [^ . 
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Table 1: Comparison between our simulation estimates of hc{T) with 
some theoretical approximations. 
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Monroe 
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WW 


BW 


WK 






(w = 0.92484) 


{uj = 0.93895) 










0.1 


3.93304(16) 


3.93307 


3.93318 


3.93069 


3.93329 


3.93330 


3.93372 


0.5 


3.6648(8) 


3.66506 


3.66561 


3.65309 


3.66611 


3.66614 


3.67589 


1.0 


3.2906(14) 


3.29303 


3.29391 


3.26843 


3.29200 


3.29261 


3.31764 


1.5 


2.7258(14) 


2.73243 


2.73396 


2.70401 


2.73094 


2.73176 


2.75099 


2.0 


1.696(2) 


1.71629 


1.71872 


1.69490 


1.71492 


1.71499 


1.71512 



TES: Tomographic entropic sampling [6] 
MHZ: Miiller-Hartmann and Zittartz [8] 
WW: Wu and Wu [l2] 
BW: Blote and Wu [B] 
WK: Wang and Kim [16] 
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Figure 9: Accessible and extrapolated configuration classes of a system of size L = 10. 
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Figure 10: Staggered magnetization per site as a function of h at T = 0.2, for L = 10 to 
30. The absolute uncertainties are plotted in the inset. 
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Fi gurc 111 Staggered susceptibility per site as a function of h 'At T = 0.02, for = 10 to 
30. The highest peaks correspond to the largest system sizes. The absolute uncertainties are 
plotted in the inset. 
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Fi gure 12: Critical field determination using finite size scaling analysis: hc{T = 0.02) 
3.98666(3). 
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Figure 13: Phase diagram of the Ising AF on the square lattice. Comparison between 
simulation and the theoretical expression of Miiller-Hartmann and Zittartz. Asterisks denote 
points obtained varying h with T fixed; circles denote points obtained varying T, with h fixed. 
The error bars of our results are smaller than the symbols. 
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